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I. INTRODUCTION 

The Gutzwiller variational wave function together with 
the Gutzwiller approximation (GA)EJ is a widely used 
approach in order to deal with Hubbard-type models. 
Originally introduced in order to explore the possibility 
of ferromagnetism within the Hubbard model (see e.g., 
Ref. H and references therein) its popularity resides in the 
fact that it captures correlation effects like the band nar- 
rowing already on the variational level. More recently 
the GA has beeiL also used for realistic band structure 
computations Since in the Hubbard model one has 
a competition between derealization, from the hopping 
of the charge carriers, and localization, from the onsite 
interaction U, the idea is to apply a projector to a given 
Slater determinant which reduces the number of doubly 
occupied sites. Within the GA one has to minimize an 
energy functional which is composed of a renormalized 
kinetic term and the interaction energy UD, where D 
denotes the concentration of doubly occupied sites. 

On the other hand, mean-field theories, like Hartree- 
Fock (HF), are usually only the first step in a many-body 
computation and it is often desirable to include the effect 
of fluctuations within the random-phase approximation 
(RPA). In case of HF this has been achieved by numer- 
ous techniques (for an overview see e.g., Ref. |^), however, 
the development of a similar scheme in the GA has been 
a long-standing problem of the condensed-matter many- 
body community. The major step in this direction was 
the reformulation of the GA by Kotliar and Ruckensteip. 
(KR) within the so-called four slave-boson approach.El 
This method maps the physical hole (or particle) into 
products of fermion and boson operators where the lat- 
ter additionally label the occupancy of the site. At the 
saddle-point level the bosons are replaced by their mean- 
field values and one recovers the GA energy functional 
showing its underlying mean-field character. 



The KR slave-boson formulation offers the possibil- 
ity of going beyond the Gutzwiller result bg the inclu- 
sion of transversal spin degrees of freedom.13 Moreover, 
in principle, it provides a controlled scheme of including 
fluctuations beyond the mean-field solution. However, 
the expansion of the KR hopping factor z SB is a highly 
nontrivial task, both with respect to the proper normal 
ordering of the bosons and with respect ta the correct 
continuum limit of the functional integral.u Expansions 
around the slave-boson saddle point have been performed 
for homogeneous systems in Refs. ||,|l0] in order to calcu- 
late correlation functions in the charge and longitudi- 
nal spin channels. Furthermore, the optical conductivity 
in the paramagnetic regime of the Hubbard model has 
been calculated in Ref. [n|. A severe difficulty in this ap- 
proach is the fact that the KR choice for .the hopping 
factor did not lead to controlled sum rules.E3 Moreover, 
to our knowledge, this approach has not been extended 
to symmetry broken states due to the complexity of the 
computation. 

Recently, two of us have presented a computation of 
RPA fluctuations on top of GA states (GA+RPA) H Our 
approach borrows, ideas from well developed techniques in 
nuclear physics Jlil and RPA fluctuations are obtained in 
the small oscillation limit of a time-dependent Gutzwiller 
approximation. Since response functions are derived for 
systems with completely unrestricted charge and spin dis- 
tributions GA+RPA is suitable also for the calculation 
of charge excitations on solutions with inhomogeneous 
textures. A key point of the GA+RPA approach is the 
proper determination of the time dependence of the varia- 
tional double occupancy parameter. We have adopted an 
antiadiabatic approximation in the sense that the double 
occupancy adjusts instantaneously to the time evolution 
of the single particle densities. In this context our ap- 
proach can be viewed as a.-generalization of the Fermi 
liquid analysis of Vollhardt.E-3 
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In this paper we use the GA+RPA to compute vari- 
ous correlation functions in the one-band Hubbard model 
and compare with exact diagonalization and HF+RPA 
results. 

This paper is organized as follows. In Sec. [n]we present 
the formalism. We concentrate on the case of the one- 
band Hubbard model although generalization to more 
complicated models is straightforward. From the ex- 
pansion of the GA energy functional up to second order 
in the densities we demonstrate how the RPA response 
functions can be calculated and demonstrate that stan- 
dard sum rules are obeyed. The method is illustrated 
in Sec. Ill for the two-site Hubbard model which can be 



treated analytically. Finally, in Sec. IV we compare the 
GA+RPA excitation spectra with exact diagonalization 
and HF+RPA results respectively. 



II. MODEL AND FORMALISM 

We consider the one-band Hubbard model 



(1) 



1J,C 



where c£1 destroys (creates) an electron with spin a at 

site i, and — c* ia Ci tC! . U is the onsite Hubbard re- 
pulsion and Uj denotes the hopping parameter between 
sites i and j. 



A. Gutzwiller approximation 



charge and spin distribution of the SD. The \ii >a act as 
local chemical potentials and are determined within the 
GA by the infinite dimension prescription that the diag- 
onal charges are not renormalized: 

(*|n<, ff |¥) = (SD\n ha \SD). 

As a—cesult one obtains for the GA energy func- 
tionalffl 

E GA [p, D}=^2 t«W^ + A, (4) 



z GA = 



y/(l - Pit + Dj)(Pii,a - Pi) + \jDj{pii-a - Pi) 
\J Pm,ct(1 — Pii,a) 



which is a functional of the density matrix p^ , 



(SPlcj^a^lSP) and the double occupancy parameters 
Pi. We denote the set of all matrix elements {pij,a} and 
{Pi} by p and P. Notice that in this paper we do not 
consider spin-canted solutions which would have density 
matrix elements (SP\Cj a Ci y(7 > \SP) + 1 for a + 1 a' . How- 
ever, transverse spin degrees of freedom can be straight- 
forwardly incorporated within the spin-rotationally in- 
variant slave boson formulation.Q 

We will denote by |\I/o) the particular wave function 
of form Eq. (^) that minimizes the energy. In order to 
obtain the associate stationary solution , one has 
to minimize E GA with respect to the double occupancy 
parameters P and the density matrix p, where the latter 
variation has to be constrained to the subspace of Slater 
determinants by imposing the projector condition p 2 = 



In its original formulation, the GA yields an approx- 
imation for the energy of a uniform paramagnetic sys- 
temlll Only in the late 80's this approach has been con- 
sistently generalized to an unrestricted Slater determi- 
nant within the Kotliar and Ruckenstein slave-boson ap- 
proach!] The same unrestricted Gutzwiller energy func- 
tional has been obtain by Gebhardli3 exploiting the fact 
that the GA becomes the exact solution of the Gutzwiller 
variational problem in the limit of infinite spatial dimen- 
sions .113 

In Gcbhard-is4ormulation the variational wave function 



is written as 



i i 



Ui = cxp ^-7i7j iiT n ia - ^ fJ-i^n-i^ 



(2) 
(3) 



where K t = (^\UiUi\^). In Eq. (|) \SP) de- 
notes a Slater determinant which already incorporates 
the Hartree contribution of the local interactions and 
which has to be determined variationally. The solution 
of the variational problem in the limit of infinite dimen- 
sions turns out to be the GA generalized for an arbitrary 



5{E<* A [p,P]-tr[A(p*-p)]} = 0, 



(5) 



where A denotes the Lagrange parameter matrix-. It is 
convenient to define a Gutzwiller Hamiltonian:li3o 



h ija [p,P] = 



dE GA 

dpjia ' 



(6) 



which is also a functional of p and P. Variation of Eq. (jfj) 
with respect to the density matrix leads to: 



h - pA - Ap + A = 0. 



(7) 



The Lagrange parameters can be eliminated^ and to- 
gether with the variation with respect to P we obtain 
the self consistent GA equations: 



[h,p] 
dE GA 
~dP~ 



0. 



0. 



(8) 
(9) 



The first equation can be solved by diagonalizing both 
the Gutzwiller Hamiltonian and the density matrix by a 
linear transformation of the single-particle orbital basis: 



(10) 
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leading to hP = <5 M „e„. Moreover, the diagonalized den- 
sity matrix p = p^ has eigenvalue 1 below the Fermi 
level and eigenvalue above it. We use Greek letters to 
denotate any state of this particular basis and the nought 
indicates evaluation in the saddle-point. Additionally we 
denotate states below the Fermi level as hole (h) states 
and the states above the Fermi level as particle states 

(P)- 

Notice that in this base (p^ ') 2 = p^ ' is trivially satis- 
fied, acts as a projector onto the hole states of the 
saddle-point Slater determinant in the space of the den- 
sity matrices whereas ct*- -* = 1 — p*- -* corresponds to the 
projector onto particle states. 

The diagonalization of Eq. (||) has to be supplemented 
by the minimization of the Gutzwiller energy with re- 
spect to the double occupancy parameters of Eq. (||) . A 
convenient method for this purpose in order to obtain in- 
homogeneous GA solutions has been discussed in Ref. |2(]. 
Notice that the unrestricted variational procedure with 
respect to charge and (or) spin degrees of freedom pre- 
vents the occurrenqe-.of the Brinkmann-Rice transition 
towards localizationcj which has already been shown in 
Ref. ^ for Neel-type antiferromagnetism. 



B. Derivation of the RPA equation 

Before starting our analysis it is convenient for later 
use to define the GA effective operator: 



GA = Y,^ A A^ + h.c), 



(11) 



•J" 



where of£ = q i3a o tj ^ and q ija = 1 if i = j and q lja 



z GA z GA otherwise. 



In order to derive the RPA equation we introduce a 
small time-dependent external field added to Eq. (Q): 



(12) 



ij a 



with fij t<T (t) = fij tCr {0)e~ lu)t . As a consequence |$), 
\SD), and the variational parameters acquire a time de- 
pendence and an additional term appears in the energy 
functional Eq. (§): 

Ef A [p,D](t) = (*(t)\H f (t)\*(t)) 

= Efc^' + M. (13) 

ijcr 

where f£j A = q^fij^. 

The time dependent field induces small amplitude os- 
cillations of D and p around the GA saddle-point: 



D = D^+SD(t), 
p = p { ' 



- »M->-6p(!t). 



(14) 
(15) 



The density and double occupancy fluctuations are con- 
strained by the following requirements: 



i) At all times p is constrained to be the one-body den- 
sity matrix associated with a Slater determinant. This 
can be achieved by imposing: 



P = P 



(16) 



ii) The double occupancy is assumed to have a much 
faster dynamics than the density matrix so that it can be 
treated antiadiabatically. As a consequence, 3D adjusts 
instantaneously to the evolution of the density matrix via 
the condition 



8E GA [ Pl D] 
dSD, 



0. 



(17) 



In fact Eq. ( |l7j ) constitutes the basic hypothesis of the 
present formalism which is necessary in order to derive 
an effective Gutzwiller interaction between particles (see 
below) . We expect this approximation to be accurate for 
sufficiently low-energy excitations. At high energies one 
can check the accuracy and the limits of validity of this 
approximation by comparing with exact diagonalization, 
as done in the following sections. Surprisingly, it turns 
out to be accurate at least up to energies of the order of 
the Mott-Hubbard gap. 

As in any small amplitude approximation, we start by 
expanding the GA energy [Eqs. (0) and (filf)] around the 
saddle-point. The first part, Eq. (|J), is needed up to 
second order in the density and double occupancy devi- 
ations: 

E GA [p,D] = E +tr{h°5p) 

+ J2^[^ lZ f A + z GA ^z GA \8p ]%a 



Here Eq denotes the saddle-point (mean-field) energy and 
the trace includes sum over spins. We have used the 
following abbreviations for the z-factor expansion 



GA 



1Z 

~dD, 
ld 2 z GA 



2 



dz 



GA 



dpu y o 

2 i \ " 



Pi; 



d 2 z GA 



2 3D 



d 2 z GA 



-y 

2 ,„ dpu.a'dpr, 



a 1 

0~Pii,a>Sp 



(19) 
(20) 



To proceed further it is convenient to cast the second 
order expression in matrix form 

E GA [p,D] = E +tr(h S P ) + ^6p ji , a Ll° H 6p lkt<T , 

+ -SDiKijSDj +SD i S l M a Sp lka , (21) 



4 



where the matrix multiplications imply the Einstein sum 
convention and the definitions for the matrices L, K and 
S follow immediately from Eqs. ©, @ and ©. The 
nonzero matrix elements are given by: 



11,(7 



^ti.jj - do - do 



J ij,ii ~ Ll 3 



dpu,a 



(22) 



A i7 — / an an \Pij,v ~T~ Pji,v) 1 T 1 3i 



v dDi dDj 



f " J d Did p. 



i i .a 



— Uj y J a „ a (Pij.cr' + Pji,cr') * 7^ Jl 



2,cr 



i ¥=3- 



Notice the formal similarity between Eq. (|T]) and an 
electron-boson problem where particle-hole excitations 
interact with a bosonic degree of freedom in the place of 
D. The matrix K plays the role of a double occupancy 
stiffness and S a double occupancy-electron interaction. 

We can integrate out the D fluctuations using the an- 
tiadiabaticity condition Eq. (jlj). First, we express 5Di 
in terms of the density fluctuations via 



SDi = -{K ^tjSjMcrSpi^a, 



(23) 



which finally yields an expansion of the energy as a func- 
tional of dp alone E[p] = E GA [p, D{p)\. 



E[p] = E Q + tr(h°Sp) 



(24) 



-Spji^lLo - SlK 1 S Q ]f J a kl 5pik^' 



Notice that we could have derived Eq. (g4|) also within 
the KR slave-boson approach. The corresponding trans- 
formations for the derivatives are given in the Appendix. 

The matrix (Lo — SqKq 1 Sq) can be considered as an 
effective interaction kernel between particle-hole excita- 
tions in the GA. For the paramagnetic regime this kernel 
reduces to the xuuasiparticle kernel of Vollhardt's Fermi 
liquid analysis.Ej Interestingly, the off-diagonal elements 
of the matrices Kij, L^ kl and Si^ia can induce inter- 
site interactions between the GA quasiparticles. This is 
in contrast with conventional HF theory of the Hubbard 
model which is purely local. 

The expansion of Ef A [p,D], Eq. (||), is needed up to 



first order only, since it is linear in the external field: 

If' 



Ef A [p,D] = F Q + tr(f GA 6p) 



EO r dqij a 
pj%aJv>°jr„ — d P kk ° 

... , ^Pkka' 
ijkaa 

En f d<lij<y 
Pii,<rJij,° 

ijko 



6D k . 



(25) 



Here Fo = f GA p^ describes the energy contribution 
when the system would be frozen at the saddle-point level 
and we used the fact that qij a does not depend on off- 
diagonal densities. As before, the double occupancy fluc- 
tuations can be eliminated through Eq. ([23]). We define 
Ef A [p] = Ef A [p,D{p)] and 



8Ef A [p] 

dpjia 



(26) 



In this paper we will restrict to density-density and 
current-current response functions with the current op- 
erator given by: 



and 



(ij) 



jij — ] tij ( c ia c jcr C^ a Ci a ). 



(27) 



When only densities are involved fij a is diagonal in the 
site index, only q ii(J = 1 is present and the last two terms 
in Eq. (|2^) vanish. If currents are involved it is easy to 
show directly from Eq. (^5|) that the last two terms also 
vanish in the absence of currents in the ground state, i.e. 

7 _ fGA 
Jo — Jo ■ 

Now, we proceed in analogy with the nuclear physics 
treatment of effective mean-field theories itussihich the 
interaction potential is density dependent £3rB Indeed 
Eq. (|24|) can be viewed as the energy expansion of an 
effective mean-field theory with the only difference that 
part of the density dependence is due to the GA hopping 
renormalization factors in the kinetic part of the Hamil- 
tonian. The advantage of this method with respect to 
other methods (e.g., equation of motion or diagrammatic 
methods) is that the present derivation is solely based on 
the knowledge of an energy functional associated with a 
Slater determinant which is precisely what the Gutzwillcr 
approximation provides. 

The density matrix of an effective meaa-field theory of 
this kind obeys the equation of mot ion rWH 

ihp=[h[p] + f(t),p}, (28) 

where we have defined an effective Gutzwiller Hamilto- 



fkjcrlp] 



dE 
dp 



(29) 



j" 7 
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which depends on densities only. At the saddle-point, we 
have ho = h[p^] = ho- The RPA is obtain by consider- 
ing the limit of small amplitude fluctuations in Eq. (28). 

It is convenient to define the four sub-sectors of the 
fluctuations of the density matrix using the projector 
properties of the density matrix discussed above: 



S P hh 


= p<®8ppl°\ 


(30) 


Sp pp 


= a^SpaW, 


(31) 


Sp hp 


= p^Spa^, 


(32) 


Sp ph 


= a^SppW. 


(33) 



The Slater determinant condition Eq. (16) implies that 
the fluctuations, Eqs. (|30|)-(|3^), are not independent. In 
fact, in terms of the fluctuations Eq. (flq) reads: 



(34) 



two particle-hole excitations. A and B are related to 
M = (L - S$Ko l S ) via 

A p h,p'h' = ( £ p — eh)8 P p'5hh' 

ij<T,nm<T f 
ij<T,mncr' 

and the transformation amplitudes ipi^ a {v) have been de- 
fined in Eq. @. 

To lowest order, we can now linearize Eq. ( pSj ) retaining 
only ph and hp matrix elements: 



Sp = p^Sp + 5pp w + {5pf 



8h ~ 
ihdp=[h Q ,Sp] + [—Sp + f,p^}, 

where we use the shorthand notation 



(41) 



Projecting Eq. ( p4| ) onto the hh and pp sector of the 
saddle point Slater determinant yields 



hh 



5p 
5p pp 



-(I + Sp hh )- 1 5p hp 5p ph m -5p hp Sp ph , (35) 
(1 - 5p pp )- 1 6p ph 6p hp w 5p ph 5p hp , (36) 



where the result on the right is valid in the small am- 
plitude limit. Thus it turns out that pp and hh density 
projections are actually quadratic in the ph and hp ma- 
trix elements. Therefore, on computing h from Eqs. ( [24] ) 
and (^9|) one should be aware of the fact that the term 
tr(h°Sp) = J^fi £ m/ 9 mm (which is first order in the pp and 
hh density projections) yields a quadratic contribution 
in the ph and hp matrix elements: 



tr(h°6p) 



E 

p 

E 

ph 



CpSppp + E e h$Phh 
h 

(fp — £h)PphPhp- 



(37) 



In addition one can neglect the pp and hh matrix element 
in the last term of Eq. (p4|). Thus, up to second order in 
the particle-hole density fluctuations, one obtains for the 
energy expansion Eq. ( |24| ) 



E[p]=E + l -{5 P hp ,8p ph ) 



A B 

B* A* 



8p hp 



(38) 



Here the so called RPA matrices A and B are given by 

dhph 



A p h,p'h' — (fp — th)&pp'5hh' 

dhp h 



dpp<h' 



B 



ph,p f h' 



dph'p' 



(39) 
(40) 



where the matrix A contains matrix elements between 
particle-hole excitations whereas the matrix B is com- 
posed of matrix elements between the ground state and 




(42) 



Then from Eqs. @, © and (|l|) one obtains the fol- 
lowing linear response equation: 



B* A* J HuJ ( -1 ) } ( Z hp ) 



fph 
fhp 



This inhomogeneous equation can be solved by inverting 
the matrix on the left-hand side which yields a linear 
relation between the external field and the change in the 
density: 



Sp = R(u)f. 



(44) 



We are now in the position to compute the response of a 
one particle observable 

O = E' : "<.^'Vt'> + /l - c -)> 

ija 

since in analogy with Eqs. (|l^) and ([[3]) its time evolution 
is given by 

(*(t)|0|*(t)) = E(Oa-W + h.c), (45) 

ija 

and the time evolution of p is known from Eq. (]44|). 

The linear response matrix R(u>) has poles at the eigen- 
frequencics of the eigenvalue problem corresponding to 
Eq. (p|) with / = 0: 



A B 

B* A* 



— Ml n 



(J -\)}(?S)=°- <«> 



Here Ml n = E n — Eq denote the excitation energies of 
the system. In analogy with the HF+RPA approxima- 
tion the vacuum of these excitations is not the old start- 
ing GA state \^o) but a new state with both Gutzwiller 



6 



type correlations and RPA ground-state correlations. We 
denote this state by |$q) an d the corresponding exited 
states by |$ n ). 

The matrix R can be written in the following Lchmann 
representation: 



R(u) P h,p'h' — 



n>0 



uj — fl n + ie ijj + £l n + it 



I p'h' 1 ph 



(47) 



In analogy with the HF+RPA method, we introduce the 
following notations: 



(oi4%» 



X 



ph> 



Ka h \n) = Yfc. 



(48) 
(49) 



The states \n) are not true excitations of the system but 
represent auxiliary notational objects. Roughly speak- 
ing, they can be thought of as RPA states without the 
Gutzwiller projector. For example |0) is the analog of the 
state \SD) but at RPA level (it contains RPA ground- 
state correlations but lacks Gutzwiller correlations). We 
will call them unprojected RPA states. The eigenvector 

(X^ , Y$ ) can be identified with the particle- hole and 
hole-particle components of the unprojected RPA excited 
state \n) with respect to the unprojected RPA ground 
state |0). 

Schematically the four states are related in the follow- 
ing way: 



\SD) 
RPA I 
10} 



l*o) 
I RPA 

l*o>, 



where P indicates Gutzwiller projection. 

Within the above formalism, it is straightforward to 
evaluate the current-current correlation function. The 
real part of the optical conductivity consists of a Drude 
part at lo — and a regular part for u> > 0: 

a(uj) = Dd(uj) +n — 0(uj - (E n - E )). 



n>0 



En — Eq 



With the above approximations and notations 
(*„Ua|*o) = (n\3° A \0), 



(50) 



where the matrix element on the right can be evalu- 
ated using Eqs. @, @, @, and @. Obviously 
the matrix elements within cr(u>) are renormalized by 
the GA hopping factors whereas R(u>) does not contain 
such renormalization. Thus the latter quantity does not 
correspond to a physical response function within the 
GA+RPA approach. 

The Drude weight D can be obtained from the f-sum 
rule (see S* -1 in next section) 



duia(io) 



~-^{T a )GA, 



(51) 



where the kinetic energy in a-direction {T^GA is evalu- 
ated within the GA. 

In practice, for computational purposes one can use 
all the standard formulas of the HF+RPA scheme by 
substituting true operators with GA effective operators 
and excitations by the unprojected excitations \n). 

The matrix elements ($o|0|<I> n ) = (0\O GA \n) can be 
used to characterize a given RPA excitation. Specific ex- 
amples which will be considered below are the transition 
density 



8n™ = (0\ni\m), 
and the transition current 

8& EE (0|jg», 



(52) 



(53) 



which can be interpreted as follows. Consider a wave 
packet 

|-0C*)> = exp(-iE t)\<P ) +r}exp(-iE m t)\$ m ), 

consisting of a small admixture rj of an exited state m 
to the ground state. For example this can be the result 
of an excitation of the mode m by an appropriate weak 
external perturbation. The time dependent expectation 
value of the charge is then given by: 



(V(t)|n*M*)) = (0\fii\Q)+vSn. 



h.c, 



and an analogous expression holds also for the current. 
Here (0|ni|0) = (SD\hi\SD) since one-particle densities 
are not renormalized by the RPA. We see that the transi- 
tion charges and currents are proportional to the ampli- 
tude of the time dependent fluctuation that would occur 
at frequeiacy Q m if the state m is excited by a weak per- 
turbationll3 (see also Ref. |3|) . 



C. Sum rules 

Sum rules form a very important tool in the theory 
of collective excitations. In many cases they allow us to 
calculate global properties in a simple way and there- 
fore they are useful in testing different approximation 
schemes. In general, a sum rule is related to the k-th 
moment of the excitation strength distribution produced 
by a single-particle operator O (see e.g., Ref. |24f ) : 



S k =J2(E n -E ) k \(^ n \O\9 ) 



(54) 



Within the present scheme we have: 

S k = J2( E n-E ) k \(n\O GA \0)\ 2 , (55) 

n 

and we restrict to current or density operators for O. The 
energy sum rule S 1 can be written as a double commu- 
tator 

S 1 = 52(E n -Eo)\{* n \0\Hfo)\ 2 = £<*o|[G, [S,O]\9 ). 

n=l 

(56) 
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In analogy with the derivation by ThoulessEl, one can 
show that the sum rule Eq. (^6|) is satisfied if the left- 
hand side is evaluated at the GA+RPA level and the 
right-hand side is calculated using the GA ground state 
wave function. The same applies for the S 1-1 sum rule, 
which in case of the optical conductivity discussed above 
corresponds to the f-sum rule. 

In the following, we consider as an example the first 
moment of the charge-charge correlation function by set- 
ting O = rii = ^2 a n ia for some lattice site R;. It is 
straightforward to evaluate the double commutator and 
we find for the corresponding sum rule 

Si = -2^2(E m - E )\(0\ ni \m)\ 2 = 2(T) GA , (57) 



where (T)ga denotes the kinetic energy evaluated within 
the GA. 

The sum rules of Eqs. (|5|) and (|5?]) provide a first en- 
couraging argument that the unrestricted GA could im- 
prove the description of charge fluctuation with respect 
to the corresponding HF method. This is based on the 
fact that the GA kinetic energy is already renormalized 
on the mean-field level. In Fig. |l| we compare the exact 
kinetic energy with unrestricted GA and HF results for 
various hole concentrations in a Hubbard model (4x4 
lattice) with nearest neighbor hopping £y = —t. The 
GA Slater determinants have been obtained using the 
method described in Ref. Notice that for this small 
system it is in general not a problem to find the true 
mean-field ground state via the variational procedure. 
We usually performed several runs starting from differ- 
ent initial configurations and checked the stability of the 
resulting states by adding some noise to the solutions. 
These are in general characterized by an inhomogeneous 
charge distribution except for the closed shell configura- 
tions. 

For small values of U/t there is almost perfect agree- 
ment between the GA method and exact results. In 
this limit, where kinetic effects dominate the correla- 
tion part, HF overestimates the value of (— T) since the 
corresponding quasiparticle hopping between sites i and 
j is described by the bare matrix element ty. On the 
other hand, for U/t — 8 the large HF onsite renormal- 
ization (corresponding to an overestimate of the spin- 
polarization) is the origin for the kinetic energy to be 
lower than the exact result. In contrast, the values of 
(— T) in the GA approximation correctly reproduce the 
exact result, especially in high-doping regime, where the 
spin density is reduced in large parts of the lattice. It 
follows that the first moment of a density-density corre- 
lation function will be more accurate in GA+RPA than 
in the HF+RPA approximation. The same holds for the 
S^-sum rule for the optical conductivity. 

To summarize this section, the idea of our method is to 
supplement the Gutzwiller approximation with REA fluc- 
tuations analogous to the HF+RPA approach. Ej Since 
the GA provides a much better initial saddle-point than 
HF one can expect that the fluctuation corrections within 
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FIG. 1: Kinetic energy per site of a 4 x 4 Hubbard model 
for U/t = 4 (a) and U/t = 8 (b) with periodic boundary 
conditions as a function of hole doping. Filled squares: exact 
result; Circles: Unrestricted HF approximation; Diamonds: 
Unrestricted GA. 

GA+RPA will allow us for a more accurate description 
of correlation effects than the HF+RPA approach. In 
the remaining sections we analyze in detail small data 
cases to test the domain of applicability of the method 
and finally show some applications in larger systems. 

III. TWO-SITE HUBBARD MODEL 

In order to demonstrate the method developed above, 
we consider in the following the two-site Hubbard model 
which can be solved exactly and can be studied analyt- 
ically in both the GA+RPA and HF+RPA approxima- 
tions. 

Exact ground-state energy and double occupancy at 
half filling (i.e., two particles) are given by 

U . 



and 



(«T n i) 



El 



2 At 2 



El" 



respectively. The exact optical conductivity displays one 
transition between the ground state and an excited state 
with energy En — U resulting in the excitation energy 



Eu-Eq 



2 (1 + ^ 



{At/Uf 



The corresponding matrix element of the current opera- 
tor is 



\mu)\ 2 = 



16< 4 



(E 2 +^) 
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Upon minimizating the GA energy functional of the half- 
filled two-site model one finds a paramagnetic solution 
below 

E^/t = 8(>/2-l)«3.31 1 
and a Neel-type state with m\ = —m\ where 

m = (n i)T ) - (th,i) 7^ 0, 

for U > U^.f t . Within HF theory the corresponding crit- 
ical value is U^ r f t /t = 2. 

Clearly the transition in cither case is non-physical 
since it does not occur in the exact solution. In this 
sense the increase of U®f t with respect to U^ r [ t is in fa- 
vor of the GA since it extends the parameter range of 
the right singlet paramagnetic solution. At large U, dis- 
regarding the non-physical broken symmetry, the Neel- 
type state in GA allows the system to reduce the double 
occupancy and at the same time prevents the occurrence 
of the Brinkmann-Rice (BR) transition towards localiza- 
tion at Ubr = 8t. 

Since the analytic expressions for the symmetry-broken 
regime become quite lengthy we restrict the derivation 
below to the paramagnetic case. In this limit the mean- 
field part of the energy is given by 



tr(ft-op) = t(l - u 2 ) }XPpp,v ~ Phh.a), 



(58) 



which defines the diagonal Gutzwiller Hamiltonian in 
Eq. (^). The hole (particle) state is the bonding (an- 
tibonding) state and u — U/(8t). 
The GA kinetic energy reads as 

E km = -2t(l - u 2 ), 

and the expansion of the GA energy functional leads to 
[see Eqs. ©] 



T (7(7 

ii.ii 



J (7(7 

it, 3 3 

J (7(7 



it (2 + 3u 2 - u 4 ) 



1 -u 2 



St 
~ u' 
8tu 2 



K, 



l-u 2 
2tu i ^ j, 
32t 



K ln = - 



l-u 2 ' 
32tu 2 



Si Ha 



l-u 2 

m 

I6tu 2 



l-u 2 
—Atu i =/= j 



i ^ 3, 



i + 3, 



(59) 
(60) 

(61) 
(62) 
(63) 

(64) 
(65) 

(66) 
(67) 



One of the peculiarities of the present approach is the 
appearance of onsite interactions for quasiparticles with 



the same spin (Lff u ^ 0). These interactions do not 
occur in the standard RPA since they would violate the 
Pauli exclusion principle but appear here because of the 
density dependence of the effective interaction between 
particles. Notice also that many of the matrix elements 
would diverge at the BR transition if it were not hidden 
by the spin-density wave (SDW) transition. 

Eliminating the double occupancy fluctuations with 
help of the antiadiabatic condition Eq. ( |l7| ) , the following 
interaction matrix is obtained: 



M aa - 



Atu 2 (3 



Ml 



M aa '- 

M aa '- 

13,13 



l-u 2 
8tu 



l-u 27 

i ± j, 

i^j, 

M aa '- = 
13,31 



-tvf 



i + 3- 



(68) 

(69) 

(70) 
(71) 
(72) 



Remarkably, intersite interactions vanish except for the 
appearance of a new interaction term between off- 
diagonal charges (M|^A Using Eqs. ( p5| ) and ( |36"| ) one 
can show that these new off-diagonal interactions do not 
contribute to the RPA matrices and the expansion of the 
energy reads as: 

2 

E = E GA + Y / [U s (Sm i ) 2 + U c (5p l ) 2 }, (73) 
i=i 

where the charge- and spin-interaction coefficients are 
given by 



ir = u< 2 -"X 1 + tt >t, 

l-u 
l + u 



(74) 
(75) 



and naturally coincide with the Landau parameters Fq , 
Fg derived by Vollhardt in Ref. [l|. 
The RPA-matrices read as 



A = 
B = 



AE + U+ U- 
U~ AE + U+ 

u+ u- 

TJ- U+ 



(76) 
(77) 



with the GA particle- hole excitation energy AE = 2t(l — 
u 2 ) and U ± = U C ±U S . 

Upon diagonalizing the RPA problem one obtains the 
eigenvalues: 



4 = AE[AE + 2{U + ± U~)], 



(78) 



which correspond to a singlet (u>+) and a magnetic ex- 
citation respectively. The former is the charge 
excitation which contributes to the optical conductiv- 
ity whereas the latter can be identified with the Gold- 
stone mode driving the transition from the paramag- 
netic to the SDW state. This transition occurs at 
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uP_ = = AE + 4U S so that the transition to the 
symmetry-broken state is only determined by the spin 
interaction U s . Notice that for the HF approximation 
we have U c = U /4 = —U s so that in this case the tran- 
sition occurs at U^.f t /t = 2 as stated above. 

When we expand the RPA charge excitation energy of 
the GA approach for small U/t we obtain 

u>\ = 2t(2t+ U) + ^U 2 +0{U A ) (GW + RPA)(79) 
u\ = 2t(2t+U) (HF + RPA), (80) 

which has to be compared with the expansion of the exact 
excitation 

u, 2 exact = 2t(2t+U) + ±U 2 + <D(U 3 ), (81) 

Thus the RPA corrections to the Gutzwiller approxima- 
tion partially includes higher order contributions in U 
which are not contained in the HF+RPA approach. 
The eigenvectors of Eq. ( |46| ) are given by 

Xf = l/2a±(l + w±/(A£7)), (82) 

Xf = ±l/2o±(l + w±/(AfO), (83) 

if = l/2a±(l-w ± /(A£)), (84) 

Yf = ±l/2a ± (l-w±/(A£)), (85) 

and the normalization factor is ot± = AE/(2w±). 

We are now able to compute the RPA double occu- 
pancy by evaluating the corresponding correlation func- 
tion as J2 m= ±(0\n-\\m) (m\n i\0) leading to 

d ga A = \ + I^E(1/c + ~1/^). (86) 

Approaching the SDW transition leads to an increase of 
spin fluctuations which leads to a divergent D RPA at the 
transition point due to the Goldstone mode w_ — ► 
(see inset of Fig. ^). From the double occupancy we 
can compute the corrections to the ground state— jea- 
ergy using the coupling constant integration trickeflEj, 
i.e. E int = Jq dxD RPA {x) which yields 

Egf +RPA = -2t + U/2 + 2t(-2+ 

+ v/2-(l- U ) 2 + v /2-(l + w) 2 ) .(87) 

We thus obtain the GA+RPA ground-state energy which 
is displayed in Fig. |^ together with the corresponding 
HF+RPA and the exact result. 

Naturally for such a small system any mean-field treat- 
ment will result in large errors, however, due to the in- 
creased value of U^ r tt and the extra contributions in U 
discussed above GA+RPA performs much better than 
HF+RPA. In general, RPA-corrections overshoot the 
exact energy which becomes significant when one ap- 
proaches the non-physical SDW transition. 
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FIG. 2: Comparison of the GA+RPA, HF+RPA and exact 
results for the ground-state energy Eo of the two-site Hubbard 
model. The inset shows the double occupancy as a function 
of U /t within the same approaches. 



For larger systems where the broken symmetry mean- 
field can correspond to long-range order or quasi long 
range order in the exact solution the agreement is much 
better as shown in Ref. [l3|. In addition the performance 
improves with increasing dimensionality as in any mean- 
field+RPA theory. 

Finally it is straightforward to check that standard 
sum rules are obeyed. For example the current opera- 
tor matrix elements between ground state and the charge 
(+) and the magnetic (— ) excitation are given by: 

l(0|i GA |-)| 2 = 0, (88) 
l(0|i GA |+)| 2 = t(l-«>+, (89) 

as a result one thus obtains that the f-sum rule 
l(0|j GA |+)| 2 A+ + (1/2)E%£ = is satisfied within the 
GA+RPA approach. 



IV. RESPONSE FUNCTIONS IN THE ID AND 
2D HUBBARD MODEL 

In this section we apply our method to the calculation 
of response functions in the ID and 2D Hubbard model. 
In the first part, we show that already on the saddle-point 
level, the GA yields rather accurate excitation energies 
as compared to HF. Inclusion of RPA corrections then 
lead to an additional redistribution of spectral weight in 
the correlation functions, which is demonstrated in the 
second part by a detailed comparison of exact diagonal- 
ization results with GA+RPA and HF+RPA. Especially 
this comparative study is intended for an a posteriori 
justification of the antiadiabaticity assumption Eq. fll7|). 
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A. Optical conductivity of the half-filled Hubbard 
model in the Gutzwiller approximation 

After Ref. [l| was published we become aware that the 
RPA residual interaction for the SDW vanishes in the 
channel relevant for the optical conductivity (zero mo- 
mentum and odd parity). This is a pathology of the 
Hubbard model at half-filling and occurs both in the 
HF+RPA and in the GA+RPA, whereas it does not oc- 
cur for more complicated models (like multiband Hub- 
bard) or symmetry-broken ground-states, like polarons 
or stripes. As a consequence the optical conductivity 
ct(lu) is the same on the GA+RPA and GA level. It 
is nevertheless quite instructive to examine <j(ui) within 
the mean-field approximation since this demonstrates the 
better starting saddle-point of GA in comparison to HF, 
with respect to the one-particle excitation energies. For 
this purpose, we compare in the following the GA and HF 
optical conductivity with numerical results and present 
the results for charge-charge correlation functions where 
the mentioned pathology does not occur and instead 
GA+RPA introduces non trivial corrections to the dy- 
namical response functions. 

Fig. H displays <t(lu) for a half- filled Hubbard chain and 
U/t = 3 convoluted with a Lorentzian L(u>) — e/[ir(e 2 + 
lu 2 )] and e = O.lt. For both GA and HF, the ground 
state is characterized by long-range SDW order and the 
regular part of (t(ui) is given by 



a(ou) 



A 2 
2lu 2 



5? \ 



j 2 - A 2 



(90) 



where t e ff contains the z-factor renormalization in case 
of the GA {t^tfj = tZfZi), t^j^ = t and 3? denotes the real 
part. The SDW gap in HF is related to the onsite mag- 
netization A HF — 2U\S Z \, whereas, within the Kotliar- 
Ruckenstein formulation of the GAB it is determined by 
the difference in the local spin-dependent Lagrange mul- 
tipliers A GA = Aj — A|. It is quite interesting to observe 
that the paset of excitations coincides rather well in the 
DDMRGE3 and GA approaches, whereas the HF gives a 
gap which is by far to large. However, although GA leads 
to excellent results for the gap energies it turns out that 
the corresponding intensity is overestimated. This has 
the consequence that most of the high-frequency evolu- 
tion of <j(lu) is compressed close to the threshold, whereas 
the DDMRG shows a much broader spectrum. We have 
checked the broadness of the exact spectrum by perform- 
ing exact diagonalization in small clusters. In fact from 
Eq. (^0|) one obtains that the large frequency tail for GA 
and HF behaves as a(ui ^> A) ~ 1/ui 3, whereas-the cor- 
rect field theoretical result is o(lu ^S> A) ~ l/wnll 

To summarize, the optical conductivity results are the 
same at mean-field or mean-field plus RPA, and, there- 
fore, no corrections are introduced by our method in this 
(rather pathological) case. As compared to the DDMRG 
results, GA performs much better than HF since it re- 
produces the onset of excitations with a better accuracy. 
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FIG. 3: Optical conductivity of the ID Hubbard chain (120 
sites) for U/t = 3. Solid: DDMRG; Dotted: GA; Dashed- 
dotted: HF. DDMRG data by courtesy of E. Jeckelmann. 
ct(ui) has been convoluted with a Lorentzian with width e = 
O.lt. 



However, the width of the spectrum is underestimated in 
both mean-field approaches. 



B. Comparison with exact results 

In order to compare the GA+RPA approach with exact 
results, we investigate the onsite density-density response 
function 



S C ( W ) = E E K*mN*0)| 2 <^ ~ (Era - Eq)). (91) 
i m>0 

In this case GA and GA+RPA solutions are different 
since the pathology discussed in the previous subsection 
is not present. 

The following sum rule is obeyed [see Eq. (|57|)1 : 



dujujSc{w) = —(T) 



GA- 



(92) 



In Fig. U we show S c (lu) for a half- filled Hubbard 
chain with 14 sites calculated with exact diagonalization, 
GA+RPA, and HF+RPA. For U/t = 3 the lowest energy 
excitation is at to ~ lAt (exact), oj ~ 1.3f (GA+RPA) 
and lu - 2. It (HF+RPA), respectively, so that GA+RPA 
is much more accurate than the standard HF+RPA ap- 
proach. Also the higher energy excitations computed 
with GA+RPA are in remarkable agreement with the 
exact results. Moreover, the oscillator strength of the 
two lowest excitations coincides rather well with the in- 
tensity obtained with exact diagonalization. The small 
high-energy features between lu ~ bt and lu ~ It present 



11 



in the exact result do not show up in the GA+RPA cor- 
relations. As a consequence, GA+RPA slightly overes- 
timates the intensity between lu ~ 3t and oj ~ 5t, since 
from (T)ca ~ (T) exact, the sum rule requires the inte- 
grated spectral weight to be approximately the same in 
both GA and the exact result. 

The accuracy of the GA+RPA approach is also re- 
markable with respect to the fact that the underlying 
mean-field solution (GA) corresponds to a SDW state, 
whereas the exact solution in ID does not show long- 
range order. However, it is well known that correlation 
functions decay slowly and thus are quasi-long ranged. 
Hence, only for very low energies one expects disagree- 
ment due to this problem and since excitations in the 
system under consideration are gapped this discrepancy 
does not really show up. 

Fig. [|(b) reveals that GA+RPA provides a better de- 
scription of the low-energy excitations than HF+RPA, 
even at larger values of U ft. In addition, we show in 
Fig. [|(b) the charge-charge correlations for the bare GA. 
In this case, the corresponding excitations are located in 
a narrow energy window on the low-energy side of the 
exact spectrum with rather incorrect oscillator strength. 
Thus RPA corrections induce the broadening and shift 
of excitations to higher energies with a simultaneous re- 
distribution of intensity. This can also be deduced from 
the inset to Fig. ||(b) which shows the evolution of the 
integrated spectrum J Q S c (y)dv for GA and GA+RPA. 
From the sum rule Eq. (^2|) it is obvious that both GA 
and GA+RPA approach the same integrated spectral in- 
tensity but with the GA+RPA spectral evolution broad- 
ened and shifted to higher energies with respect to the 
bare GA. 

In systems with large dimensionality we expect our 
approach to improve, since the GA is an exact solution 
of the Gutzwiller variational problem in infinite dimen- 
sions. Furthermore, quite generally mean-field theories 
become better as the dimensionality is increased. Fig. |5| 
displays S c (uj) for a half-filled 4x4 system (i.e. 16 parti- 
cles) with U/t = 10. The lowest prominent energy peak 
in the exact solution occurs at ujf£ in = 8.4t and a second 
bunch of excitations starts at U/t w lOt . . . lit decaying 
in intensity towards higher energy. The lowest energy ex- 
citation within GA+RPA (w%f+ RPA = 8.7t) is remark- 
ably close to the exact value in contrast to HF+RPA 
where the lowest peak appears at u)^J n = 9.8t. More- 
over, the center of high energy excitations in the exact 
solution is represented by two peaks in the GA+RPA 
spectrum at 9.7t and 11. It whereas they are shifted to 
slightly higher energy within HF+RPA. It is interesting 
to observe that also in this energy range the GA+RPA 
method gives a better (although rather crude) approx- 
imation than HF+RPA despite the expected failure of 
the antiadiabaticity condition Eq. ( |l7| ) for energies larger 
than the Mott-Hubbard gap. 

We finally would like to demonstrate the importance 
of an accurate mean-field solution for the quality of the 
RPA excitation spectrum. For this purpose Fig. o, shows 




FIG. 4: Charge correlation function S c (u}) for a half-filled 
Hubbard chain (14 site) in case of U/t = 3 (a) and U/t — 6 
(b). For clarity the curves for HF+RPA and GA (in b) have 
been shifted upwards. The inset in (b) shows the integrated 
weight as a function of frequency for GA and GA+RPA, re- 
spectively. 



2 - 



3_ 



exact 




U/t=10 


GA+RPA 






HF+RPA 


il 
i! 


I N p= 16 


i\ 


\ i 1 




1 V.' 


\i 


V \ A 














; 8 


11 


14 



03/1 



FIG. 5: Charge correlation function Sc(uj) for the half-filled 
4x4 Hubbard model (16 particles) in case of U/t — 10. 




FIG. 6: Charge correlation function S c (w) for the 4x4 Hub- 
bard model and 10 particles in case of U/t — 4 (a) and 
U/t= 10(b). 



S c (ui) for a 4 x 4 system and 10 particles. In the case 
of U/t = 4 [Fig. | (a)] HF+RPA and GA+RPA give a 
good approximation to the lowest excitation, both with 
respect to the oscillator strength and the energy. For 
U/t = 4 and 10 particles (corresponding to a closed shell 
configuration) the underlying mean-field solution are ho- 
mogeneous with respect to the charge and do not show 
SDW order. However, whereas the homogeneous GA so- 
lution remains a stable saddle-point also for large values 
of U/t, the HF solution becomes instable with respect to 
a disordered charge and spin texture. This obviously has 
dramatic consequences for the dynamical properties as 
shown in Fig. ^|(b). In fact, the GA+RPA spectrum still 
shows a remarkable agreement with the exact solution in 
contrast to the RPA on top of the inhomogeneous HF 
solution. 

Notice that the GA+RPA oscillator strength for the 
charge-charge correlations is distributed in a much better 
way than for the optical conductivity (see Fig. ||). To 
some extent, this can be attributed to the constraints 
imposed by the sum rules. Both correlation functions 
satisfy a sum rule which involves the kinetic energy on the 
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right hand side [Eqs. (|l]) and J57j)], however, in case of 
the optical conductivity the high-energy states contribute 
much less to the sum rule than in the case of S c (u>) due 
to the w factor. Therefore, the high-energy part of S c (oj) 
is much more constrained to be accurate than the high- 
energy part of a. This argument assumes that at least the 
excitation energies are approximately correct, otherwise 
one can satisfy the sum rule by a compensation between 
the incorrect excitation energies and the incorrect matrix 
elements as it occurs within the bare GA for the charge 
correlation function [Fig. ^(b)] . 

Finally, it is remarkable that the GA+RPA dynamical 
correlation functions perform rather well even at energies 
much larger than the charge gap (Figs. ^J, |5| and ||) where 
the antiadiabatic condition Eq. ( |17| ) is not expected to 
hold. To some extent this may be due to the constraints 
imposed by sum rules. It shows that at least for some cor- 
relation functions (charge-charge) the domain of applica- 
bility of our theory may be much wider than expected. 



V. CONCLUSION 

In summary, we have presented a method for the cal- 
culation of dynamical correlation functions in the Hub- 
bard model, based on the Gutzwiller approximation. 
Our method obeys well-behaved sum rules and we have 
demonstrated that it is suitable for practical compu- 
tations: excitation energies compare remarkably bet- 
ter with exact diagonalization results than the related 
HF+RPA approach. Moreover, since the performance of 
any RPA computation crucially depends on the quality 
of the underlying mean-field solution we conclude from 
our analysis that the GA provides a much better starting 
point for this purpose than HF. 

The dynamical matrix has been calculated as a 
quadratic expansion of the GA energy functional in the 
densities, which allows us to construct the RPA eigen- 
vectors and eigenvalues. Essential assumption for carry- 
ing out this expansion is the antiadiabaticity condition 
Eq. (0). Despite the fact that charge and double occu- 
pancy dynamics seem to be governed by different time 
scales it would be desirable to relax the antiadiabatic ap- 
proximation and to treat the double occupancy dynamics 
explicitly. This is in principle possible via the Kotliar- 
Ruckenstein slave-boson scheme, however, attempts in 
this direction render difficult due to the hopping factor 
expansion. 

Compared to numerical methodsc3c3 our approach can 
be pushed to much .larger systems. Our experience on 
modeling real dataEIH is that often finite size effects are 
more severe than the inaccuracies introduced by mean- 
field+RPA approaches. A more recent approach for dy- 
namical properties consists of mapping the problem onto 
quantum impurity models (dynamical mean-field the- 
ory) which becomes exact in the limit of large dimen- 
sions.E3 This has enormously increased our understand- 
ing of these systems. However, on making the limit of 
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large dimensions important parts of the physics are lost. 
For example all acoustic like collective behavior, like spin 
waves disappears. On the other hand, these collective ef- 
fects are naturally captured in our approach. 

The GA+RPA formalism can be also applied to multi- 
band Hubbard models, which are relevant for a more 
quantitative analysis of excitations in the cuprate high- 
T c superconductors. It has been shown recently that the 
GA provides an excellent starting point to describe the 
physics of stripes in cuprates including the behavior of 
inconmensurability, chemical potential and some trans- 
port experiments with doping.E3 In this context it is very 
important to make an RPA analysis on top of GA states 
since within the HF approximation one obtains a ground 
state which does not correspond to experiment. Work in 
this direction is in progress. 
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APPENDIX A: RELATIONSHIP TO 
KOTLIAR-RUCKENSTEIN SLAVE BOSON 
APPROACH 

In the Kotliar^jRuckenstein slave-boson approach to the 
Hubbard model] the original Hilbert space is enlarged by 
introducing four subsidiary boson fields e*, Sjj, s^i, and 
di for each site R^. These operators stand for the an- 
nihilation of empty, singly occupied states with spin up 
or down, and doubly occupied sites, respectively. Since 
there are only four possible states per site, these boson 
projection operators must satisfy the completeness con- 
straints: 



and 



e]e; 



d]di = 1, 



(Al) 



d]d, 



(A2) 



In the saddle-point approximation, all bosonic operators 
are treated as numbers and the resulting effective one- 
particle Hamiltonian H KR describes the dynamics of par- 
ticles where the hopping amplitude between states (i,cr) 
and (j,cr) is renormalized by a factor zf B zf B with 



-SB 



{£ + s l-a) V2 {e%Si t<T + Si-vdi) (d 2 + s? >0 .) 
The total energy of H KR is given by 



-1/2 

(A3) 



E SB = 



(A4) 



t z SB z SB o ■ 

ij.a i 

which has to be minimized (i) with res pect t o th e bosonic 
fields within the constraints Eqs. (Al) and ([A4 ) and (ii) 
with respect to p within the subspace of Slater determi- 
nants. 

The slave-boson energy functional E SB is a function 
of 4N boson variables where N is the number of lattice 



sites. Sin ce these bosons obey the constraints Eqs. (Al) 
and (A2) one can eliminate 3N of them which leads to 
the Gutzwiller energy E GA when one keeps the double 
occupancy variable D = d 2 for each lattice site. Thus 
the expansions of E SB and E GA are connected via the 
transformation 
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2 
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de 2 



(A5) 



and the derivatives have to be taken at the saddle-point. 
Upon inserting this transformation in Eqs. ([H]) and ( |20| ) 
leads to an analogous energy expansion than Eq. ( p4| ) but 
now within the KR scheme. For paramagnetic solutions 
this corresponds to the analysis done in Ref. nO. 
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